######################################################
# CreatePovertyVariables.R
# This script creates school-specific poverty estimates 
# using WorldPop data. Note that the code for the raster merge
# cannot be executed without coordinates (which have been removed to comply with IRB)
# The raster merge code has been included below for reference however. 
# note this file requires the working directory to be set to "./Merge and Clean Data/"
# Contact Ryan Jablonski, r.s.jablonski@lse.ac.uk with questions
#
# Log
# Created 2018
# Edited for APSR replication 17 August 2023 by Ryan Jablonski
######################################################


rm(list=ls(all=TRUE))
library(groundhog)

groundhog.library(maps,'2023-04-01')

library(foreign)
groundhog.library(sp,'2023-04-01')
library(maptools)
gpclibPermit()
library(classInt)
groundhog.library(rgeos,'2023-04-01')
library(dismo)
groundhog.library(rgdal,'2023-04-01')
library(raster)

#p4s = CRS("+proj=longlat +datum=WGS84")
#r.proj=CRS("+proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")

#r=raster(".\\input\\mwi11povcons125.tif")
#r.unc=raster(".\\input\\mwi11povcons125uncert.tif")


schooldata.lc=read.csv("./output/Schools.forLC.withpop.csv", stringsAsFactors=FALSE)
schooldata.mp=read.csv("./output/Schools.forMP.withpop.csv", stringsAsFactors=FALSE)
# 
# school.points <- SpatialPointsDataFrame(coords = schooldata.lc[, c("school_longitude","school_latitude")], data = schooldata.lc,
#                                         proj4string = p4s)
# school.points <- spTransform(school.points, r.proj)
# 
# data.extract=extract(r, school.points)
# data.extract.unc=extract(r.unc, school.points)
# 
# r.df=(as.data.frame(r))
# r.unc.df=(as.data.frame(r.unc))
# 
# r.df$upper=r.df$mwi11povcons125+r.unc.df$mwi11povcons125uncert
# r.df$lower=r.df$mwi11povcons125-r.unc.df$mwi11povcons125uncert
# 
# hist(r.df$mwi11povcons125)
# hist(r.df$lower)
# hist(r.df$upper)
# 
# 
# 
# data.malawi.merge = cbind(as.data.frame(school.points), data.frame(poverty_proportion=data.extract), data.frame(poverty_proportion_unc=data.extract.unc))
# write.csv(data.malawi.merge[,c("school_id", "poverty_proportion_unc", "poverty_proportion")], "./input/worldpop_poverty_school_lc.csv")

poverty.lc=read.csv("./input/worldpop_poverty_school_lc.csv")
schooldata.lc$poverty_proportion=poverty.lc[match(schooldata.lc$school_id, poverty.lc$school_id),"poverty_proportion"]
schooldata.lc$poverty_proportion_unc=poverty.lc[match(schooldata.lc$school_id, poverty.lc$school_id),"poverty_proportion_unc"]

write.csv(schooldata.lc, "./output/Schools.forLC.withpoppoverty.csv")

#MP
# school.points <- SpatialPointsDataFrame(coords = schooldata.mp[, c("school_longitude","school_latitude")], data = schooldata.mp,
#                                         proj4string = p4s)
# school.points <- spTransform(school.points, r.proj)
# 
# data.extract=extract(r, school.points)
# data.extract.unc=extract(r.unc, school.points)
# 
# r.df=(as.data.frame(r))
# 
# data.malawi.merge = cbind(as.data.frame(school.points), data.frame(poverty_proportion=data.extract), data.frame(poverty_proportion_unc=data.extract.unc))
# write.csv(data.malawi.merge[,c("school_id", "poverty_proportion_unc", "poverty_proportion")], "./input/worldpop_poverty_school_mp.csv")

poverty.mp=read.csv("./input/worldpop_poverty_school_mp.csv")
schooldata.mp$poverty_proportion=poverty.mp[match(schooldata.mp$school_id, poverty.mp$school_id),"poverty_proportion"]
schooldata.mp$poverty_proportion_unc=poverty.mp[match(schooldata.mp$school_id, poverty.mp$school_id),"poverty_proportion_unc"]

write.csv(schooldata.mp, "./output/Schools.forMP.withpoppoverty.csv")

